#include "kernel/cell.hpp"
//#include "yocto/math/fcn/zfind.hpp"

void Cell:: steadyState()
{
	Solution       &Stmp    = solutions[ iS_tmp ];
	Solution       &dS      = solutions[ idS ];
	const Solution &S_in    = solutions[ iS_in ];
	const Solution &S_out   = solutions[ iS_out ];
	
	std::cerr << "-------- Steady State --------" << std::endl;
	//--------------------------------------------------------------------------
	//-- user Em
	//--------------------------------------------------------------------------
	std::cerr << "-- Em=" << Em*1000.0 << " mV" << std::endl;
	std::cerr << "-- Capacity=" << Capa << std::endl;
	
	//--------------------------------------------------------------------------
	//-- compute passive leaks
	//--------------------------------------------------------------------------
	library.reset_permeabilities();
	const Real zeta_m = Em * F_over_RT;
	dS.zero();
	dS.addLeak(zeta_m, Volume, S_in, S_out );
	
	
	//-- check valid K+
	if( dS[ "K+" ] >= 0 )
		throw exception("Invalid [K+] values!");
	
	//-- check valid chloride
	if( dS[ "Cl-"] >= 0)
		throw exception("Invalid [Cl-] values!");
	
	//--------------------------------------------------------------------------
	//-- sum everyone but K+
	//--------------------------------------------------------------------------
	Real ZLeaks = 0;
	for( size_t i=dS.size; i>0; --i )
	{
		const Species &species = dS(i);
		if( species.name != "K+" )
		{
			ZLeaks += species.Z * dS[i];
		}
	}
	std::cerr << "ZLeaks=" << ZLeaks << std::endl;
	if( ZLeaks <= 0 )
		throw exception("Invalid passive flux can't match Em");
	
	const Real lambda_K_manip = dS["K+"]; assert( lambda_K_manip < 0 );
	const Real lambda_K_match = - (2.0*ZLeaks)/3.0;
	std::cerr << "lambda_K_manip = " << lambda_K_manip << std::endl;
	std::cerr << "lambda_K_match = " << lambda_K_match << std::endl;
	const Real KPf      = lambda_K_match / lambda_K_manip;
	std::cerr << "Correcting by " << KPf << std::endl; 
	library["K+"].Pf = KPf;
	
	//--------------------------------------------------------------------------
	//-- compute matching enzymatic activities: recompute passive dS
	//--------------------------------------------------------------------------
	dS.zero();
	dS.addLeak(zeta_m, Volume, S_in, S_out );
	
	//--------------------------------------------------------------------------
	//-- compute NaK coefficient
	//--------------------------------------------------------------------------
	pNaK->V = 1;
	const Real rho_NaK1 = pNaK->addRate(Stmp,0.0,zeta_m,S_in,S_out,*this);
	const Real lambda_K = dS[ "K+" ];
	std::cerr << "-- rho_NaK1=" << rho_NaK1 << std::endl;
	std::cerr << "-- lambda_K=" << lambda_K << std::endl;
	pNaK->V = -0.5 * (lambda_K/rho_NaK1);
	if( pNaK <= 0 )
		throw exception("Invalid Potassium Concentrations!");
	
	//--------------------------------------------------------------------------
	//-- compute AE->V
	//--------------------------------------------------------------------------
	pAE->V = 1;
	const Real lambda_Cl = dS["Cl-"];
	const Real rho_AE1   = pAE->addRate(Stmp,0.0,zeta_m,S_in,S_out,*this);
	std::cerr << "-- lambda_Cl=" << lambda_Cl << std::endl;
	std::cerr << "-- rho_AE1  =" << rho_AE1 << std::endl;
	const Real Vmax = -lambda_Cl;
	pAE->V = (Vmax/rho_AE1);
	if( pAE->V <= 0 )
		throw exception("Invalid Chloride Concentrations!");
	
	//--------------------------------------------------------------------------
	//-- compute NHE->V since rho_NHE = rho_AE
	//--------------------------------------------------------------------------
	pNHE->V = 1;
	const Real rho_NHE1 = pNHE->addRate(Stmp,0.0,zeta_m,S_in,S_out,*this);
	std::cerr << "-- rho_NHE1= " << rho_NHE1 << std::endl;

	pNHE->V = (Vmax/rho_NHE1);
	std::cerr << "-- V_AE  = " << pAE->V << std::endl;
	std::cerr << "-- V_NHE = " << pNHE->V << std::endl;
    std::cerr << "-- Ratio = " << (pAE->V)/(pNHE->V) << std::endl;
    
	computeRate( dS, 0.0, Em, S_in);
	std::cerr << "dS0=" << std::endl;
	dS.display(std::cerr);
	std::cerr << "--------" << std::endl << std::endl;
	
}
